home *** CD-ROM | disk | FTP | other *** search
/ QRZ! Ham Radio 8 / QRZ Ham Radio Callsign Database - Volume 8.iso / mac / files / ant_nec / nec81tar.z / nec81tar / intrp.f < prev    next >
Text File  |  1991-05-13  |  12KB  |  443 lines
  1. C $TITLE: 'INTRP'
  2. C $NOFLOATCALLS
  3. C
  4.       SUBROUTINE INTRP (X,Y,F1,F2,F3,F4)
  5. C
  6. C     INTRP USES BIVARIATE CUBIC INTERPOLATION TO OBTAIN THE VALUES OF
  7. C     4 FUNCTIONS AT THE POINT (X,Y).
  8. C
  9.       REAL*8 X,Y,DX,DY,XS,YS,XZ,YZ,XX,YY
  10.       COMPLEX*16 F1,F2,F3,F4
  11.       COMPLEX*16 A,B,C,D,FX1,FX2,FX3,FX4,P1,P2,P3,P4,
  12.      1A11,A12,A13,A14,A21,A22,A23,A24,A31,A32,A33,A34,A41,A42,A43,A44,
  13.      2B11,B12,B13,B14,B21,B22,B23,B24,B31,B32,B33,B34,B41,B42,B43,B44,
  14.      3C11,C12,C13,C14,C21,C22,C23,C24,C31,C32,C33,C34,C41,C42,C43,C44,
  15.      4D11,D12,D13,D14,D21,D22,D23,D24,D31,D32,D33,D34,D41,D42,D43,D44
  16. C**
  17.       REAL*4 DXA,DYA,XSA,YSA
  18.       COMPLEX*8 AR1,AR2,AR3,EPSCF,ARL1,ARL2,ARL3
  19. C**
  20.       COMMON/GGRID/AR1(11,10,4),AR2(17,5,4),AR3(9,8,4),EPSCF,DXA(3),
  21.      1 DYA(3),XSA(3),YSA(3),NXA(3),NYA(3)
  22.       DIMENSION NDA(3), NDPA(3)
  23.       DIMENSION A(4,4),B(4,4),C(4,4),D(4,4),ARL1(1),ARL2(1),ARL3(1)
  24.       EQUIVALENCE(A(1,1),A11),(A(1,2),A12),(A(1,3),A13),(A(1,4),A14)
  25.       EQUIVALENCE(A(2,1),A21),(A(2,2),A22),(A(2,3),A23),(A(2,4),A24)
  26.       EQUIVALENCE(A(3,1),A31),(A(3,2),A32),(A(3,3),A33),(A(3,4),A34)
  27.       EQUIVALENCE(A(4,1),A41),(A(4,2),A42),(A(4,3),A43),(A(4,4),A44)
  28.       EQUIVALENCE(B(1,1),B11),(B(1,2),B12),(B(1,3),B13),(B(1,4),B14)
  29.       EQUIVALENCE(B(2,1),B21),(B(2,2),B22),(B(2,3),B23),(B(2,4),B24)
  30.       EQUIVALENCE(B(3,1),B31),(B(3,2),B32),(B(3,3),B33),(B(3,4),B34)
  31.       EQUIVALENCE(B(4,1),B41),(B(4,2),B42),(B(4,3),B43),(B(4,4),B44)
  32.       EQUIVALENCE(C(1,1),C11),(C(1,2),C12),(C(1,3),C13),(C(1,4),C14)
  33.       EQUIVALENCE(C(2,1),C21),(C(2,2),C22),(C(2,3),C23),(C(2,4),C24)
  34.       EQUIVALENCE(C(3,1),C31),(C(3,2),C32),(C(3,3),C33),(C(3,4),C34)
  35.       EQUIVALENCE(C(4,1),C41),(C(4,2),C42),(C(4,3),C43),(C(4,4),C44)
  36.       EQUIVALENCE(D(1,1),D11),(D(1,2),D12),(D(1,3),D13),(D(1,4),D14)
  37.       EQUIVALENCE(D(2,1),D21),(D(2,2),D22),(D(2,3),D23),(D(2,4),D24)
  38.       EQUIVALENCE(D(3,1),D31),(D(3,2),D32),(D(3,3),D33),(D(3,4),D34)
  39.       EQUIVALENCE(D(4,1),D41),(D(4,2),D42),(D(4,3),D43),(D(4,4),D44)
  40.       EQUIVALENCE(ARL1,AR1),(ARL2,AR2),(ARL3,AR3),(XS2,XSA(2)),
  41.      1 (YS3,YSA(3))
  42.       DATA IXS,IYS,IGRS/-10,-10,-10/,DX,DY,XS,YS/1.,1.,0.,0./
  43.       DATA NDA/11,17,9/,NDPA/110,85,72/,IXEG,IYEG/0,0/
  44.       IF (X.LT.XS.OR.Y.LT.YS) GO TO 1
  45.       IX=INT((X-XS)/DX)+1
  46.       IY=INT((Y-YS)/DY)+1
  47. C
  48. C     IF POINT LIES IN SAME 4 BY 4 POINT REGION AS PREVIOUS POINT, OLD
  49. C     VALUES ARE REUSED
  50. C
  51.       IF (IX.LT.IXEG.OR.IY.LT.IYEG) GO TO 1
  52.       IF (IABS(IX-IXS).LT.2.AND.IABS(IY-IYS).LT.2) GO TO 12
  53. C
  54. C     DETERMINE CORRECT GRID AND GRID REGION
  55. C
  56. 1     IF (X.GT.XS2) GO TO 2
  57.       IGR=1
  58.       GO TO 3
  59. 2     IGR=2
  60.       IF (Y.GT.YS3) IGR=3
  61. 3     IF (IGR.EQ.IGRS) GO TO 4
  62.       IGRS=IGR
  63.       DX=DXA(IGRS)
  64.       DY=DYA(IGRS)
  65.       XS=XSA(IGRS)
  66.       YS=YSA(IGRS)
  67.       NXM2=NXA(IGRS)-2
  68.       NYM2=NYA(IGRS)-2
  69.       NXMS=((NXM2+1)/3)*3+1
  70.       NYMS=((NYM2+1)/3)*3+1
  71.       ND=NDA(IGRS)
  72.       NDP=NDPA(IGRS)
  73.       IX=INT((X-XS)/DX)+1
  74.       IY=INT((Y-YS)/DY)+1
  75. 4     IXS=((IX-1)/3)*3+2
  76.       IF (IXS.LT.2) IXS=2
  77.       IXEG=-10000
  78.       IF (IXS.LE.NXM2) GO TO 5
  79.       IXS=NXM2
  80.       IXEG=NXMS
  81. 5     IYS=((IY-1)/3)*3+2
  82.       IF (IYS.LT.2) IYS=2
  83.       IYEG=-10000
  84.       IF (IYS.LE.NYM2) GO TO 6
  85.       IYS=NYM2
  86.       IYEG=NYMS
  87. C
  88. C     COMPUTE COEFFICIENTS OF 4 CUBIC POLYNOMIALS IN X FOR THE 4 GRID
  89. C     VALUES OF Y FOR EACH OF THE 4 FUNCTIONS
  90. C
  91. 6     IADZ=IXS+(IYS-3)*ND-NDP
  92.       DO 11 K=1,4
  93.       IADZ=IADZ+NDP
  94.       IADD=IADZ
  95.       DO 11 I=1,4
  96.       IADD=IADD+ND
  97.       GO TO (7,8,9), IGRS
  98. C     P1=AR1(IXS-1,IYS-2+I,K)
  99. 7     P1=ARL1(IADD-1)
  100.       P2=ARL1(IADD)
  101.       P3=ARL1(IADD+1)
  102.       P4=ARL1(IADD+2)
  103.       GO TO 10
  104. 8     P1=ARL2(IADD-1)
  105.       P2=ARL2(IADD)
  106.       P3=ARL2(IADD+1)
  107.       P4=ARL2(IADD+2)
  108.       GO TO 10
  109. 9     P1=ARL3(IADD-1)
  110.       P2=ARL3(IADD)
  111.       P3=ARL3(IADD+1)
  112.       P4=ARL3(IADD+2)
  113. 10    A(I,K)=(P4-P1+3.*(P2-P3))*.1666666667D0
  114.       B(I,K)=(P1-2.*P2+P3)*.5
  115.       C(I,K)=P3-(2.*P1+3.*P2+P4)*.1666666667D0
  116. 11    D(I,K)=P2
  117.       XZ=(IXS-1)*DX+XS
  118.       YZ=(IYS-1)*DY+YS
  119. C
  120. C     EVALUATE POLYMOMIALS IN X AND THEN USE CUBIC INTERPOLATION IN Y
  121. C     FOR EACH OF THE 4 FUNCTIONS.
  122. C
  123. 12    XX=(X-XZ)/DX
  124.       YY=(Y-YZ)/DY
  125.       FX1=((A11*XX+B11)*XX+C11)*XX+D11
  126.       FX2=((A21*XX+B21)*XX+C21)*XX+D21
  127.       FX3=((A31*XX+B31)*XX+C31)*XX+D31
  128.       FX4=((A41*XX+B41)*XX+C41)*XX+D41
  129.       P1=FX4-FX1+3.*(FX2-FX3)
  130.       P2=3.*(FX1-2.*FX2+FX3)
  131.       P3=6.*FX3-2.*FX1-3.*FX2-FX4
  132.       F1=((P1*YY+P2)*YY+P3)*YY*.1666666667D0+FX2
  133.       FX1=((A12*XX+B12)*XX+C12)*XX+D12
  134.       FX2=((A22*XX+B22)*XX+C22)*XX+D22
  135.       FX3=((A32*XX+B32)*XX+C32)*XX+D32
  136.       FX4=((A42*XX+B42)*XX+C42)*XX+D42
  137.       P1=FX4-FX1+3.*(FX2-FX3)
  138.       P2=3.*(FX1-2.*FX2+FX3)
  139.       P3=6.*FX3-2.*FX1-3.*FX2-FX4
  140.       F2=((P1*YY+P2)*YY+P3)*YY*.1666666667D0+FX2
  141.       FX1=((A13*XX+B13)*XX+C13)*XX+D13
  142.       FX2=((A23*XX+B23)*XX+C23)*XX+D23
  143.       FX3=((A33*XX+B33)*XX+C33)*XX+D33
  144.       FX4=((A43*XX+B43)*XX+C43)*XX+D43
  145.       P1=FX4-FX1+3.*(FX2-FX3)
  146.       P2=3.*(FX1-2.*FX2+FX3)
  147.       P3=6.*FX3-2.*FX1-3.*FX2-FX4
  148.       F3=((P1*YY+P2)*YY+P3)*YY*.1666666667D0+FX2
  149.       FX1=((A14*XX+B14)*XX+C14)*XX+D14
  150.       FX2=((A24*XX+B24)*XX+C24)*XX+D24
  151.       FX3=((A34*XX+B34)*XX+C34)*XX+D34
  152.       FX4=((A44*XX+B44)*XX+C44)*XX+D44
  153.       P1=FX4-FX1+3.*(FX2-FX3)
  154.       P2=3.*(FX1-2.*FX2+FX3)
  155.       P3=6.*FX3-2.*FX1-3.*FX2-FX4
  156.       F4=((P1*YY+P2)*YY+P3)*YY*.1666666667D0+FX2
  157.       RETURN
  158.       END
  159. C
  160. C
  161. C
  162.       SUBROUTINE HSFLD (XI,YI,ZI,AI)
  163. C     HSFLD COMPUTES THE H FIELD FOR CONSTANT, SINE, AND COSINE
  164. C     CURRENT ON A SEGMENT INCLUDING GROUND EFFECTS.
  165.       REAL*8 XI,YI,ZI,RHOSPC
  166.       COMPLEX*16 EXK,EYK,EZK,EXS,EYS,EZS,EXC,EYC,EZC,ZRATI,ZRATI2,
  167.      1 T1,FRATI,RRV,RRH,ZRATX,QX,QY,QZ,HPK,HPS,HPC
  168.       INTEGER*4 IND1,IND2
  169.       COMMON/DATAJ/S,B,XJ,YJ,ZJ,CABJ,SABJ,SALPJ,EXK,EYK,EZK,EXS,EYS,
  170.      1 EZS,EXC,EYC,EZC,RKH,IND1,IND2,IPGND,IEXK
  171.       COMMON/GND/ZRATI,ZRATI2,FRATI,CL,CH,SCRWL,SCRWR,NRADL,KSYMP,
  172.      1 IFAR,IPERF,T1,T2
  173.       DATA ETA/376.73/
  174.       XIJ=XI-XJ
  175.       YIJ=YI-YJ
  176.       RFL=-1.
  177.       DO 7 IP=1,KSYMP
  178.       RFL=-RFL
  179.       SALPR=SALPJ*RFL
  180.       ZIJ=ZI-RFL*ZJ
  181.       ZP=XIJ*CABJ+YIJ*SABJ+ZIJ*SALPR
  182.       RHOX=XIJ-CABJ*ZP
  183.       RHOY=YIJ-SABJ*ZP
  184.       RHOZ=ZIJ-SALPR*ZP
  185. C      RH=DSQRT(RHOX*RHOX+RHOY*RHOY+RHOZ*RHOZ+AI*AI)
  186.       RH=SQRT(RHOX*RHOX+RHOY*RHOY+RHOZ*RHOZ+AI*AI)
  187.       IF (RH.GT.1.E-10) GO TO 1
  188.       EXK=0.
  189.       EYK=0.
  190.       EZK=0.
  191.       EXS=0.
  192.       EYS=0.
  193.       EZS=0.
  194.       EXC=0.
  195.       EYC=0.
  196.       EZC=0.
  197.       GO TO 7
  198. 1     RHOX=RHOX/RH
  199.       RHOY=RHOY/RH
  200.       RHOZ=RHOZ/RH
  201.       PHX=SABJ*RHOZ-SALPR*RHOY
  202.       PHY=SALPR*RHOX-CABJ*RHOZ
  203.       PHZ=CABJ*RHOY-SABJ*RHOX
  204.       CALL HSFLX (S,RH,ZP,HPK,HPS,HPC)
  205.       IF (IP.NE.2) GO TO 6
  206.       IF (IPERF.EQ.1) GO TO 5
  207.       ZRATX=ZRATI
  208. C      RMAG=DSQRT(ZP*ZP+RH*RH)
  209.       RMAG=SQRT(ZP*ZP+RH*RH)
  210. C      XYMAG=DSQRT(XIJ*XIJ+YIJ*YIJ)
  211.       XYMAG=SQRT(XIJ*XIJ+YIJ*YIJ)
  212. C
  213. C     SET PARAMETERS FOR RADIAL WIRE GROUND SCREEN.
  214. C
  215.       IF (NRADL.EQ.0) GO TO 2
  216.       XSPEC=(XI*ZJ+ZI*XJ)/(ZI+ZJ)
  217.       YSPEC=(YI*ZJ+ZI*YJ)/(ZI+ZJ)
  218. C      RHOSPC=DSQRT(XSPEC*XSPEC+YSPEC*YSPEC+T2*T2)
  219.       RHOSPC=SQRT(XSPEC*XSPEC+YSPEC*YSPEC+T2*T2)
  220.       IF (RHOSPC.GT.SCRWL) GO TO 2
  221.       RRV=T1*RHOSPC*DLOG(RHOSPC/T2)
  222.       ZRATX=(RRV*ZRATI)/(ETA*ZRATI+RRV)
  223. 2     IF (XYMAG.GT.1.E-6) GO TO 3
  224. C
  225. C     CALCULATION OF REFLECTION COEFFICIENTS WHEN GROUND IS SPECIFIED.
  226. C
  227.       PX=0.
  228.       PY=0.
  229.       CTH=1.
  230.       RRV=(1.,0.)
  231.       GO TO 4
  232. 3     PX=-YIJ/XYMAG
  233.       PY=XIJ/XYMAG
  234.       CTH=ZIJ/RMAG
  235.       RRV=CDSQRT(1.-ZRATX*ZRATX*(1.-CTH*CTH))
  236. 4     RRH=ZRATX*CTH
  237.       RRH=-(RRH-RRV)/(RRH+RRV)
  238.       RRV=ZRATX*RRV
  239.       RRV=(CTH-RRV)/(CTH+RRV)
  240.       QY=(PHX*PX+PHY*PY)*(RRV-RRH)
  241.       QX=QY*PX+PHX*RRH
  242.       QY=QY*PY+PHY*RRH
  243.       QZ=PHZ*RRH
  244.       EXK=EXK-HPK*QX
  245.       EYK=EYK-HPK*QY
  246.       EZK=EZK-HPK*QZ
  247.       EXS=EXS-HPS*QX
  248.       EYS=EYS-HPS*QY
  249.       EZS=EZS-HPS*QZ
  250.       EXC=EXC-HPC*QX
  251.       EYC=EYC-HPC*QY
  252.       EZC=EZC-HPC*QZ
  253.       GO TO 7
  254. 5     EXK=EXK-HPK*PHX
  255.       EYK=EYK-HPK*PHY
  256.       EZK=EZK-HPK*PHZ
  257.       EXS=EXS-HPS*PHX
  258.       EYS=EYS-HPS*PHY
  259.       EZS=EZS-HPS*PHZ
  260.       EXC=EXC-HPC*PHX
  261.       EYC=EYC-HPC*PHY
  262.       EZC=EZC-HPC*PHZ
  263.       GO TO 7
  264. 6     EXK=HPK*PHX
  265.       EYK=HPK*PHY
  266.       EZK=HPK*PHZ
  267.       EXS=HPS*PHX
  268.       EYS=HPS*PHY
  269.       EZS=HPS*PHZ
  270.       EXC=HPC*PHX
  271.       EYC=HPC*PHY
  272.       EZC=HPC*PHZ
  273. 7     CONTINUE
  274.       RETURN
  275.       END
  276. C
  277. C
  278. C
  279.       SUBROUTINE HSFLX (S,RH,ZPX,HPK,HPS,HPC)
  280. C     CALCULATES H FIELD OF SINE COSINE, AND CONSTANT CURRENT OF SEGMENT
  281.       REAL*8 TP,FJKX,PI8,SDK,CDK,HKR,HKI,ZP,DH,Z1,Z2,RHZ,RH2
  282.       COMPLEX*16 FJK,HPK,HPS,HPC,EKR1,EKR2,T1,T2,CONS
  283.       COMPLEX FJ
  284.       DIMENSION FJX(2), FJKX(2)
  285.       EQUIVALENCE (FJ,FJX), (FJK,FJKX)
  286.       DATA TP/6.283185308D0/,FJX/0.,1./,FJKX/0.,-6.283185308D0/
  287.       DATA PI8/25.13274123D0/
  288.       IF (RH.LT.1.E-10) GO TO 6
  289.       IF (ZPX.LT.0.) GO TO 1
  290.       ZP=ZPX
  291.       HSS=1.
  292.       GO TO 2
  293. 1     ZP=-ZPX
  294.       HSS=-1.
  295. 2     DH=.5*S
  296.       Z1=ZP+DH
  297.       Z2=ZP-DH
  298.       IF (Z2.LT.1.E-7) GO TO 3
  299.       RHZ=RH/Z2
  300.       GO TO 4
  301. 3     RHZ=1.
  302. 4     DK=TP*DH
  303. C      CDK=DCOS(DK)
  304.       CDK=COS(DK)
  305. C      SDK=DSIN(DK)
  306.       SDK=SIN(DK)
  307.       CALL HFK (-DK,DK,RH*TP,ZP*TP,HKR,HKI)
  308.       HPK=DCMPLX(HKR,HKI)
  309.       IF (RHZ.LT.1.E-3) GO TO 5
  310.       RH2=RH*RH
  311.       R1=DSQRT(RH2+Z1*Z1)
  312.       R2=DSQRT(RH2+Z2*Z2)
  313.       EKR1=CDEXP(FJK*R1)
  314.       EKR2=CDEXP(FJK*R2)
  315.       T1=Z1*EKR1/R1
  316.       T2=Z2*EKR2/R2
  317.       HPS=(CDK*(EKR2-EKR1)-FJ*SDK*(T2+T1))*HSS
  318.       HPC=-SDK*(EKR2+EKR1)-FJ*CDK*(T2-T1)
  319.       CONS=-FJ/(2.*TP*RH)
  320.       HPS=CONS*HPS
  321.       HPC=CONS*HPC
  322.       RETURN
  323. 5     EKR1=DCMPLX(CDK,SDK)/(Z2*Z2)
  324.       EKR2=DCMPLX(CDK,-SDK)/(Z1*Z1)
  325.       T1=TP*(1./Z1-1./Z2)
  326.       T2=CDEXP(FJK*ZP)*RH/PI8
  327.       HPS=T2*(T1+(EKR1+EKR2)*SDK)*HSS
  328.       HPC=T2*(-FJ*T1+(EKR1-EKR2)*CDK)
  329.       RETURN
  330. 6     HPS=(0.,0.)
  331.       HPC=(0.,0.)
  332.       HPK=(0.,0.)
  333.       RETURN
  334.       END
  335. C
  336. C
  337. C
  338.       SUBROUTINE HFK (EL1,EL2,RHK,ZPKX,SGR,SGI)
  339. C     HFK COMPUTES THE H FIELD OF A UNIFORM CURRENT FILAMENT BY
  340. C     NUMERICAL INTEGRATION
  341.       COMMON/TMH/ ZPK,RHKS
  342.       INTEGER*4 NM,NS
  343.       REAL*8 T01R,T10R,TE1R,T01I,T10I,TE1I,T11R,T20R,TE2R,T11I,
  344.      1 T20I,TE2I,T00R,T00I,RHK,ZPKX,SGR,SGI
  345.       DATA NX,NM,NTS,RX/1,65536,4,1.E-4/
  346.       ZPK=ZPKX
  347.       RHKS=RHK*RHK
  348.       Z=EL1
  349.       ZE=EL2
  350.       S=ZE-Z
  351.       EP=S/(10.*NM)
  352.       ZEND=ZE-EP
  353.       SGR=0.0
  354.       SGI=0.0
  355.       NS=NX
  356.       NT=0
  357.       CALL GH (Z,G1R,G1I)
  358. 1     DZ=S/NS
  359.       ZP=Z+DZ
  360.       IF (ZP-ZE) 3,3,2
  361. 2     DZ=ZE-Z
  362.       IF (ABS(DZ)-EP) 17,17,3
  363. 3     DZOT=DZ*.5
  364.       ZP=Z+DZOT
  365.       CALL GH (ZP,G3R,G3I)
  366.       ZP=Z+DZ
  367.       CALL GH (ZP,G5R,G5I)
  368. 4     T00R=(G1R+G5R)*DZOT
  369.       T00I=(G1I+G5I)*DZOT
  370.       T01R=(T00R+DZ*G3R)*0.5
  371.       T01I=(T00I+DZ*G3I)*0.5
  372.       T10R=(4.0*T01R-T00R)/3.0
  373.       T10I=(4.0*T01I-T00I)/3.0
  374.       CALL TEST (T01R,T10R,TE1R,T01I,T10I,TE1I,0.D0)
  375.       IF (TE1I-RX) 5,5,6
  376. 5     IF (TE1R-RX) 8,8,6
  377. 6     ZP=Z+DZ*0.25
  378.       CALL GH (ZP,G2R,G2I)
  379.       ZP=Z+DZ*0.75
  380.       CALL GH (ZP,G4R,G4I)
  381.       T02R=(T01R+DZOT*(G2R+G4R))*0.5
  382.       T02I=(T01I+DZOT*(G2I+G4I))*0.5
  383.       T11R=(4.0*T02R-T01R)/3.0
  384.       T11I=(4.0*T02I-T01I)/3.0
  385.       T20R=(16.0*T11R-T10R)/15.0
  386.       T20I=(16.0*T11I-T10I)/15.0
  387.       CALL TEST (T11R,T20R,TE2R,T11I,T20I,TE2I,0.D0)
  388.       IF (TE2I-RX) 7,7,14
  389. 7     IF (TE2R-RX) 9,9,14
  390. 8     SGR=SGR+T10R
  391.       SGI=SGI+T10I
  392.       NT=NT+2
  393.       GO TO 10
  394. 9     SGR=SGR+T20R
  395.       SGI=SGI+T20I
  396.       NT=NT+1
  397. 10    Z=Z+DZ
  398.       IF (Z-ZEND) 11,17,17
  399. 11    G1R=G5R
  400.       G1I=G5I
  401.       IF (NT-NTS) 1,12,12
  402. 12    IF (NS-NX) 1,1,13
  403. 13    NS=NS/2
  404.       NT=1
  405.       GO TO 1
  406. 14    NT=0
  407.       IF (NS-NM) 16,15,15
  408. 15    WRITE(*,18)  Z
  409.       GO TO 9
  410. 16    NS=NS*2
  411.       DZ=S/NS
  412.       DZOT=DZ*0.5
  413.       G5R=G3R
  414.       G5I=G3I
  415.       G3R=G2R
  416.       G3I=G2I
  417.       GO TO 4
  418. 17    CONTINUE
  419.       SGR=SGR*RHK*.5
  420.       SGI=SGI*RHK*.5
  421.       RETURN
  422. C
  423. 18    FORMAT (24H STEP SIZE LIMITED AT Z=,F10.5)
  424.       END
  425. C
  426. C
  427. C
  428.       SUBROUTINE GH (ZK,HR,HI)
  429. C     INTEGRAND FOR H FIELD OF A WIRE
  430.       REAL*8 CKR,SKR,RS,R,RR2,RR3
  431.       COMMON/TMH/ ZPK,RHKS
  432.       RS=ZK-ZPK
  433.       RS=RHKS+RS*RS
  434.       R=DSQRT(RS)
  435.       CKR=DCOS(R)
  436.       SKR=DSIN(R)
  437.       RR2=1./RS
  438.       RR3=RR2/R
  439.       HR=SKR*RR2+CKR*RR3
  440.       HI=CKR*RR2-SKR*RR3
  441.       RETURN
  442.       END
  443.